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In our companion paper, the physiological functions of pancreatic 3 cells were analyzed with a new (3-cell model 
by time-based integration of a set of differential equations that describe individual reaction steps or functional 
components based on experimental studies. In this study, we calculate steady-state solutions of these differential 
equations to obtain the limit cycles (LCs) as well as the equilibrium points (EPs) to make all of the time derivatives 
equal to zero. The sequential transitions from quiescence to burst-interburst oscillations and then to continuous 
firing with an increasing glucose concentration were defined objectively by the EPs or LCs for the whole set of equa- 
tions. We also demonstrated that membrane excitability changed between the extremes of a single action potential 
mode and a stable firing mode during one cycle of bursting rhythm. Membrane excitability was determined by the 
EPs or LCs of the membrane subsystem, with the slow variables fixed at each time point. Details of the mode 
changes were expressed as functions of slowly changing variables, such as intracellular [ATP], [Ca 2+ ], and [Na + ]. 
In conclusion, using our model, we could suggest quantitatively the mutual interactions among multiple membrane 
and cytosolic factors occurring in pancreatic (3 cells. 



INTRODUCTION 

In our companion paper (see Cha et al. in this issue) , 
we constructed a detailed model of a pancreatic (3 cell 
based on published electrophysiological measurements 
of ion channels or exchangers. The model successfully 
reconstructed three representative electrical activities 
in response to a varying glucose concentration ([G]): the 
quiescent states of the membrane potential (V m ) , burst- 
ing activity with alternate burst-interburst events, and 
continuous firing of action potentials. It was suggested 
that the burst-interburst cycle is generated by the inter- 
actions of channels or transporters with intracellular ions 
and/or metabolic intermediates. By applying lead po- 
tential (Vl) analysis (Cha et al., 2009), we could quantify 
contributions of individual membrane currents to the 
slow changes in V m during the interburst period, and sug- 
gested distinct ionic mechanisms underlying the burst- 
ing rhythm at different [G]. 

In our companion paper (Cha et al., 2011), the dynamic 
behavior of the model was calculated by time-based 
integration of ordinary differential equations. For ex- 
ample, [G] -dependent activity in a (3 cell was examined 
by integrating 18 differential equations until a constant 
pattern was observed at each [G]. However, this time- 
based simulation will never be more than an approxi- 
mation because of the uncertainty that always remains 
in defining the steady states even after a long integra- 
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tion time, and in discriminating different patterns of 
membrane excitability explicitly. These problems may 
be overcome if steady-state solutions of the differential 
equations are obtained with respect to continuous vari- 
ation in [G] or slowly varying cytosolic substances. To 
achieve this aim, we have applied bifurcation analysis to 
the comprehensive (3-cell model developed in our com- 
panion paper. 

Based on the steady-state solutions, the cellular re- 
sponses in a (3 cell will be explicitly described in terms 
of "mode changes" of system behavior. We focused on 
three main questions: what kind of modes contribute to 
membrane excitability in (3 cells; when does each mode 
change occur during the burst-interburst rhythm; and, 
finally, which are the important intracellular factors un- 
derlying the mode changes? Collectively, with the results 
of time-based simulations and V L analysis in our com- 
panion paper (Cha et al., 2011), the results of the bifur- 
cation analysis clarified the physiological roles of several 
intracellular factors in promoting modal changes in 
(3-cell function. The results will be compared with those 
of the bifurcation analysis to a series of simple (3-cell 
models reported by Chay, Keizer, and colleagues in the 
past few decades. 
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MATERIALS AND METHODS 

The new p-cell model as a nondrift system 
with unique solutions 

The structure and individual components of our P-cell model 
were fully described in our companion paper (Cha et al., 2011). 
In brief, the model is composed of 18 variables: V m , seven gating 
variables of ion channels, three state variables of iNaCa. four ionic 
concentrations in the cytosol ([Na + ]j, [K + ]i, and [Ca 2+ ]j) or in ER 
([Ca 2+ ] ER ), and three metabolic substrate concentrations ([ATP], 
[MgADP] , and [Re] ) . Time-dependent changes of these variables 
are described by ordinary differential equations (refer to the sup- 
plemental material in Cha et al., 2011). 

The new P-cell model is a nondrift system; that is, the full 
18-variable system has steady-state solutions to make all the time 
derivatives in the model zero simultaneously. A necessary condi- 
tion for the existence of these steady-state solutions is that all of 
the cation fluxes (Na + , K + , and Ca 2+ ) through ion channels and 
exchangers should be included in calculating the time derivatives 
of both the membrane potential (dV m / dt) and the intracellular 
ion concentrations (d[X + ]/dt) . No intracellular ion concentra- 
tion was fixed arbitrarily in the model. Additionally, to avoid un- 
limited changes of the metabolic compounds, the sum of NADH 
and NAD was set to be constant, like that of ATP and AD P. Using 
these approaches, the simulated responses of our P-cell model 
were all completely reversible. 

To obtain the unique solutions in the full system (Fig. 1), the 
redundancy in calculating four ion concentrations ([Na + ];, [Ca 2+ ]j, 
[K + ]j, and [Ca 2+ ] E R) and V m was avoided by applying the charge con- 
servation law: 



[Ca 2+ ] ER = j _ik_ (Vm - V m (0) ) - ( [Na + ], - [Na + ], (0) ) (A8) 

2 vol ER volj F 

-( [K + ], - [K + ],(0) ) - -( [Ca 2+ ] ; - [Ca 2+ ],(0) ) | + [Ca 2+ ] ER (0) , 

as derived in the Appendix. This procedure removed d [Ca 2+ ] ER / dt 
from the model, leaving 17 independent variables. The initial 
value of each variable (for example, V m (0)) are presented in 
Table SI in our companion paper (Cha et al., 2011). 

Bifurcation analysis 

In experimental studies using real 3 cells, a steady-state response 
to a new experimental condition is observed after a certain delay, 
for example, after applying a new level of [G] . The time-based sim- 
ulation is straightforward to reconstruct this experimental proto- 
col by repeating the integration of model differential equations 
with a tiny time step. In contrast to this time-based simulation, the 
bifurcation analysis directly solves multiple differential equations, 
which are responsible for the kinetic behavior of the model. 
Namely, the steady-state solutions of the model are directly ob- 
tained by setting all of the differential equations to zero and solv- 
ing them simultaneously. The solutions are usually represented in 
a bifurcation diagram, which shows changes in equilibrium points 
(EPs) and limit cycles (LCs) when one parameter is systematically 
varied on the x axis, for example, [G] in Fig. 1 or Pcav in Fig. 2. An 
EP corresponds to a steady-state point at which the system perma- 
nently stays unless any perturbations are applied, and an LC is a 
steady-state periodic solution, for example, a sustained oscillation 
in V m and substrate concentrations. In this paper, LCs are repre- 
sented with a pair of lines that indicate the maximum and mini- 
mum values during an oscillation. The stability of an EP can be 
explicitly determined by an eigenvalue of ajacobian matrix. A system 
eventually converges to a stable EP, such as the resting membrane 



potential and accompanying stable intracellular substrate con- 
centrations. A system can also stay at a hypothetical unstable EP, 
but any perturbation will cause the system to leave from the EP, 
like a ball perfectly balanced on the peak of a hill. Similarly, an 
oscillation (LC) can also be stable or unstable. 

Bifurcation refers to a behavioral change in a system, such as 
the change in the cellular electrical activity from the resting po- 
tential to the spontaneous action potential generation when the 
extracellular [G] is increased. A bifurcation is generally indicated 
with an emergence or disappearance of EPs or LCs, or with a change 
in their stability. The following kinds of bifurcation were observed in 
this study: limit point (LP) bifurcation: two EPs or LCs approach 
and annihilate each other; Hopf bifurcation: an EP loses stability 
and a new LC appears; period doubling bifurcation: the system 
switches to a new behavior with twice the period of the original 
system; and Torus bifurcation: an LC becomes unstable, and the 
system oscillates around the unstable LC. Mathematical details of 
classifying bifurcations could be skipped. In this study, the steady- 
state solutions were numerically obtained using AUTO implemented 
in XPPAUT (Ermentrout, 2002), a specific computational tool for 
bifurcation analysis. For a more extensive explanation about the 
bifurcation analysis to the cellular models, see Fall et al. (2002). 

Fast and slow decomposition of the system to determine 
membrane excitability 

Our companion paper (Cha et al., 2011) demonstrated that the 
ionic mechanism for membrane excitation is continuously modi- 
fied by the time-dependent variations in cytosolic substrates during 
burst-interburst rhythm. If the rate of these time-dependent changes 
in cytosolic substrate concentrations is slow enough to scarcely 
affect the configuration of individual action potentials, but could 
determine the burst-interburst cycle, the membrane excitability 
at a given time point might be examined by fixing the substrate 
concentrations. We applied bifurcation analysis by fixing the in- 
tracellular substrate concentrations ([S]iS) and [Ca 2+ ] ER at a specific 
time point of the burst-interburst rhythm to determine the modes 
of the membrane excitability. After fixing these concentrations, 
however, we occasionally found that a slow change in the mem- 
brane excitability still occurred solely because of the time-dependent 
changes of the ultraslow inactivation gate of I Ca v (f us ) (Fig. S2 A). 
When f us was fixed in addition to the substrate concentrations, the 
membrane excitability stayed in the same state as observed at the 
fixing moment. That is, a steady rhythm of repetitive action po- 
tentials or a steady resting potential was established depending on 
the fixed time (Fig. S2 B). Thus, fixing these variables largely 
satisfied the requirement for fast and slow decomposition to de- 
fine the membrane excitability at a given time point. Based on the 
range of time constants during bursting rhythm (Fig. S2 C), the 
inactivation state of I Na ca (L) was also classified as a slow variable 
because it has a comparable time constant to the f^. [Ca 2+ ]i oscil- 
lated rapidly over the time span of a single action potential (fast 
Ca 2+ ripple) superimposed on a slow drift of the Ca 2+ plateau dur- 
ing the burst (inset of Fig. 4 A or Fig. 2 in Cha et al., 2011). In this 
study, [Ca 2+ ] ; was treated as a slow variable because it took ^4 s for 
[Ca 2+ ]j to reach a new steady state when other slow variables in- 
cluding [Ca 2+ ] ER were fixed. As a consequence, to determine the 
mode of membrane excitability, we calculated steady-state solu- 
tions for the membrane system consisting of nine fast variables 
(V m , dcav, U CaV , rgDn qKD„ m KC a ( BK), h KC a(BK), Ei_ tota i, and L), after fixing 
nine slow variables ([ATP], [MgADP], [Re], [Na + ] i; [KTh, [Ca 2+ ]i, 
and [Ca 2+ ] ER , in addition to f^ and L) at a given time of observation. 
Refer to Fig. SI or the supplemental material in our companion 
paper (Cha et al., 2011) for definitions of individual variables. 

Online supplemental material 

Fig. SI shows the interactions between model variables in the cyto- 
solic and membrane systems. Fig. S2 shows model behaviors when 



40 Changes in membrane excitability in (3 cells 



all the intracellular concentrations ( [Sj] s) are fixed, or when Lj, was 
also fixed together with [SJs. Fig. S2 C demonstrates the time 
constants of the model variables during bursting. Fig. S3 shows a 
simulation that is started from an unstable EP at 16 mM [G]. The 
source code of this model is provided as a PDF file for XPPAUT 
(Ermentrout, 2002). The supplemental material is available at 
http://www.jgp.org/cgi/content/full/jgp.201110612/DCl. 

RESULTS 

Mode changes in cellular activity induced by varying [G] 

Extracellular glucose affects [ATP] and [MgADP] 
through ATP production pathways; the equilibrium val- 
ues of the remaining variables are readjusted through 
reciprocal interactions among the metabolic compounds, 
membrane excitation, and intracellular ion concentra- 
tions (Fig. SI). Accordingly, different patterns (or modes) 
of cellular activity are established by varying [G], as 
shown by the time-based simulation in our companion 
paper (Cha et al., 2011). In the bifurcation analysis, 
underlying mode changes are revealed explicitly by 
the steady-state solutions (EPs and LCs) of the differen- 
tial equations for all 17 variables in the model. Fig. 1 
demonstrates changes in V m , [ATP], [Ca 2+ ];, and [Na + ]; 
of EPs or LCs when [G] is varied. At [G] < 6.9 mM, one 
stable EP exists (Fig. 1, black lines), predicting that the 
cellular state will always return exactly to this EP regard- 
less of the initial condition. Namely, for [G] < 6.9 mM, 
EP defines the resting membrane potential and the cor- 
responding steady-state values of substrate concentra- 
tions. For 6.9 < [G] < 30 mM, EP becomes unstable (Fig. 1, 
red lines) indicating that spontaneous bursts of action 
potentials start to take place. Around the unstable EP, 
an LC is present represented by the yellow and blue 
lines (Fig. 1 ) , which give the peak and the most negative 
potentials of individual action potential in the case of 
stable LCs. For 6.9 < [G] < 18.84 mM, the LC is unstable 
and the action potential configuration slowly varies dur- 
ing the intermittent spike burst period. The LC converts 



from unstable to stable at [G] = 18.84 mM, predicting 
the continuous spike burst in the time-based simulation, 
(refer to Fig. 2 of Cha et al., 2011). Thus, the sequential 
transitions from quiescence to burst-interburst oscilla- 
tions and then to continuous firing with increasing [G] 
were defined objectively by the EPs or LCs obtained by 
solving the whole set of equations. 

The unstable EP for 6.9 < [G] < 30 mM (Fig. 1, red 
lines) indicates a hypothetical condition, under which 
the membrane potential and the intracellular composi- 
tion of ions and substrates can remain constant unless any 
perturbation is applied. The unstable EP (Fig. 1, red lines) 
is usually difficult to observe in experiments using real 
cells, or in time-based simulations, but its existence can 
be convinced by starting the simulation from that par- 
ticular set of values of all variables defined by the EP in the 
bifurcation analysis. An example is shown in Fig. S3 A, 
where V m remained near the initial values for ^1 s be- 
fore spontaneous activity started. This is typical behav- 
ior for any unstable EP. The "spontaneous" deviation 
from the EP is most probably a result of numerical errors 
in integration. An alternative way to put the system on 
the unstable EP in real experiments or in the time-based 
simulation might be provided by the voltage-clamp ex- 
periments. Fig. S3 B shows a protocol of clamping V m at 
the equilibrium potential given by the unstable EP for 
[G] = 16 mM. Under the voltage clamp, values of the re- 
maining 16 variables, including [ATP], [Ca 2+ ];, and [Na + ];, 
spontaneously relaxed to their equilibrium levels defined 
by the unstable EP (Fig. S3 C) . This response is expected 
for a stable EP but not for an unstable EP. We ran the 
bifurcation analysis again and found that EP became 
stable when V m , the pivotal factor in the membrane sys- 
tem, was fixed (not depicted) . When the membrane was 
released from the voltage clamp (Fig. S3 B), the EP be- 
came unstable again, and V m escaped after a similar latency 
as in Fig. S3 A. Obviously, voltage-clamp experiments 
are awaited to confirm this prediction in real cells. 
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Figure 1. Changes in EPs and LCs in the 
whole (3-cell model by varying [G] . Four bifur- 
cation diagrams showing continuous changes 
in EPs or LCs for V m (A), [ATP] (B), [Ca 2+ ]i 
(C), and [Na*]i (D) with respect to [G]. Sta- 
ble EPs, unstable EPs, stable LCs, and unsta- 
ble LCs are indicated by black, red, blue, and 
yellow lines, respectively. For LCs, the maxi- 
mum and minimum values in oscillations 
were plotted. The amplitude of oscillation 
in [ATP] and [Na + ]i was small, and the maxi- 
mum and minimum curves of LCs fused with 
each other. AUTO failed to find unstable LCs 
at [G] < 9 mM. Open circles indicate bifur- 
cation points; HB, Hopf bifurcation from a 
stable EP to an unstable EP (at 6.9 mM [G]); 
TR, Torus bifurcation from a stable LC to an 
unstable LC (at 18.84 mM [G]). 
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The [G]-EP relationships (Fig. 1, black and red lines) 
give important clues on the principal mechanisms of the 
glucose-induced signal transductions, free from the cyclic 
oscillation in time-based simulations. With increasing [G] , 
the accelerated ATP production raises the equilibrium 
level of [ATP] in a dose-dependent manner (Fig. 1 B). 
Then, the Ikatp conductance continuously decreases, and 
the resulting membrane depolarization (Fig. 1 A) in- 
creases the equilibrium level of [Ca 2+ ] ; via a fractional ac- 
tivation of IcaV: even though I S oc is deactivating (Fig. 1 C) . 
Elevated [Ca 2+ ]i accelerates ATP consumption, resulting 
in a small inflection in the equilibrium [ATP]-[G] curve 
(Fig. 1 B) . The biphasic relationship between [Na + ] i and 
[G] is determined by two factors (Fig. 1 D). The initial 
falling phase in [Na + ]; is a result of activation of the 
Na + /K + pump through an increase in [ATP] , and the ris- 
ing phase is attributed to the increase in Na + influx via 
the forward mode of Na + /Ca 2+ exchange, which is accel- 
erated by increased [Ca 2+ ];. With the saturation of ATP 
production >20 mM [G], all relationships become flat. 

Mode changes in membrane excitability induced by slow 
changes in the intracellular substrates 

Time-base simulations, as well as experimental studies, 
have led to the hypothesis that membrane excitability is 
cyclically modulated by slow changes in intracellular 
substrates during the burst-interburst rhythm. To prove 
this hypothesis mathematically, we separated the mem- 
brane system from the cytosolic factors by adopting the 
same strategy of the fast and slow decomposition of the 
system variables, as has been established in previous 
bifurcation analyses of simple (3-cell models (see Materials 
and methods). The nine slow variables were fixed at 
each time point of the bursting activity when solving the 
differential equations of the nine fast variables to obtain 
EPs or LCs in the fast membrane system. Based on the 
EPs and LCs, we can define the membrane excitability, 
as described below. 



Definition of modes of membrane excitability. In discrimi- 
nating modes of membrane excitability in our new 
(3-cell model, we found it convenient to calculate EPs and 
LCs by varying the amplitude factor (Pc a v) of Icav, a domi- 
nant current in generating action potentials. The bifur- 
cation diagram showed a typical S-shaped curve of EP in 
the Pcav-V m plane (Fig. 2 A) . The black and red lines indi- 
cate stable and unstable EPs, as in Fig. 1, and LCs are 
shown by blue (stable) or yellow (unstable) lines. The 
black vertical line in Fig. 2 A indicates the standard 
Pcav (48.9 pA mM -1 ), and the intersections with the bi- 
furcation diagram correspond to the EPs or LCs in 
the control system. 

Six distinct modes of membrane excitability (Fig. 2 A, 
top of the panel, from A to F) were defined by finding 
the values of Pcav for each bifurcation point (Fig. 2 A, gray 
lines) . Each mode has a different number of stable or 
unstable EPs and LCs, and shows specific electrophysio- 
logical characteristics in response to current injections 
into the cell (Fig. 3 and summarized in Table I). Modes 
A and F with extremely small or large Pcav are nonexcit- 
able. Mode B has two unstable EPs and one stable EP, and 
a single action potential is induced by an electrical stimu- 
lus. Mode C has two stable states with spontaneous action 
potentials (stable LC) and a resting potential (stable EP). 
Accordingly, the stable firing can be switched on and off 
by applying a depolarizing or a hyperpolarizing current 
pulse, respectively. Mode D with a stable LC shows a stable 
oscillation regardless of stimuli. Mode E is also bistable with 
stable oscillations and a depolarized resting potential. 

In the Pcav-V m plane in Fig. 2 A, the membrane excit- 
ability at the standard Pc a v is determined to be at mode 
C under the given set of slow variables. The correspond- 
ing bistable membrane excitation is obvious from the 
N-shaped I-V diagram with three intersections, drawn 
with the same set of slow variables (Fig. 2 B). The zero cur- 
rent potentials, EPi, EP 2 , and EP 3 on the I-V curve, 
correspond to the intersections of the EP curve with the 



Figure 2. Mode changes of membrane 
excitability by varying P Ca v- (A) Bifurcation 
diagrams showing continuous changes 
in V m of EPs or LCs as a function of Pc a v 
Stable EPs, unstable EPs, stable LCs, and 
unstable LCs are indicated by black, red, 
blue, and yellow lines, respectively. The 
black line for a stable EP corresponds to 
the resting membrane potential, and the 
two blue lines for a stable LC correspond 
to the amplitude of the action potentials. 
The slow variables were fixed to the follow- 
ing values: [ATP] = 2.64 mM, [MgADP] = 
0.0591 mM, [Re] = 0.61 mM, [Na + ] t = 
5.87 mM, [K + ]i = 127 mM, [Ca 2+ ]i = 0.124 pM, 

[Ca 2+ ] ER = 0.0247 mM, f^ = 0.843, and L = 0.152. Open circles indicate bifurcation points; LP, LP bifurcation; HB, Hopf bifurcation; PD, 
period doubling bifurcation. Black vertical lines with the blue italic numeral indicate control values of P Ca v (48.9 pA mM -1 ) in the (3-cell 
model. EPi, EP 2 , and EP 3 (gray dots) are the intersections of the EP curve with the black vertical line. Gray vertical lines passing through 
the corresponding bifurcation points separate individual modes, as indicated at the top. (B) Steady-state I-V relationship. Zero current 
potentials correspond to EPj, EP 2 , and EP 3 in A. 
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black vertical lines at Pc a v of 48.9 pA mM -1 in Fig. 2 A. 
The resting V m is at EP 1; and action potentials, if trig- 
gered by an appropriate stimulus, oscillate around EP 3 . 
The I-V curve, however, fails to determine if the action 
potential is repetitive or singular because of lack of in- 
formation about LCs. 

Time-dependent mode changes at 8 mM [G]. To disclose 
specific modes of membrane excitability during the 
burst-interburst rhythm at 8 mM [G], bifurcation dia- 
grams (V m vs. Pcav) were obtained using the values of [S] ;s, 
[Ch 2+ ]er, f us , and L at successive time points, measured 
from the time-based simulation shown in Fig. 4 A (gray 
line). Fig. 4 B shows bifurcation diagrams at six repre- 
sentative time points during spontaneous bursting activ- 
ity. EPs or LCs were measured at the standard value of 
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Figure 3. Representative membrane responses to current in- 
jections in different modes. The amplitudes of injected current 
pulses are indicated inside panels. Duration of the pulse was 5 ms, 
except 15 ms in Mode E. In each panel, a different value of P Ca v 
was used for simulation (units in pA mM -1 ) : Mode A, P Ca v was 30; 
Mode B, 45; Mode C, 48.9; Mode D, 60; Mode E, 65; Mode F, 69. 
The same values were used for the slow variables as in Fig. 2. The 
time axes were identical for all panels. 



Pcav (48.9 pA mM - , black vertical lines) in individual 
bifurcation diagrams and were plotted together with 
the result of time-based simulation (Fig. 4 A) . Finally, we 
determined the mode of membrane excitability at each 
time point on the basis of the definitions summarized in 
Table I and showed sequential mode changes at the top 
of Fig. 4. Mode C (bistable mode) observed at the begin- 
ning of the record is followed by mode D (stable firing 
mode) at 7.6 s, and action potentials are initiated after a 
delay. The burst was terminated by extinction of the sta- 
ble LC through a mode change from D to B (single action 
potential mode) via the temporal appearance of mode C. 
After cessation of the burst, the membrane returns again 
to mode C. It should be noted that [Ca 2+ ]; fluctuated 
rapidly during the burst, as shown in the inset in Fig. 4 A. 
Because of these rapid changes in [Ca 2+ ]i, the position 
of the EP or LC curves fluctuated slightly along the ab- 
scissa in Fig. 4 B, and we obtained bifurcation diagrams at 
the extremes of [Ca 2+ ];. The mode changes from D to C 
or from C to B occurred slighdy earlier at the local mini- 
mums of the Ca 2+ transient than at the maximums. 

Fig. 4 B indicates that a sequence of bifurcation oc- 
curred during one cycle of bursting rhythm. At 0 s, V m 
remains stable at EP^ During the initial 7.6 s, the dia- 
gram shifts leftward with reference to the standard Pcav, 
and EPi and EP 2 approach each other. At 7.6 s, EPi co- 
alesces with EP 2 and disappears (LP bifurcation of EPs) , 
and thus the membrane system shifts to the stable LC, 
corresponding with spontaneous action potentials. It takes 
about another 3 s for V m to move from EPi to LC because 
of the very small inward current charging the membrane 
capacitance. Once spontaneous oscillations are initi- 
ated, the diagram shifts rightward until the stable LC 
disappears by an LP bifurcation at 15.5 s. Finally, V m re- 
turns to the stable EP 1; and the whole cycle of events re- 
peats again. Here, it should be noted that the movement 
of the bifurcation diagram is entirely the result of time- 
dependent changes in the slow intracellular factors. 

Which slow variables are responsible for termination 
of the burst? 

Termination of the burst occurs from the disappear- 
ance of the stable LC (LP bifurcation) , accompanied by 

TABLE I 

Membrane excitability 



Mode 


EPs or LCs 


Physiology 


A 


1 sEP 


Nonexcitable mode 


B 


1 sEP, 2 uEP 


Single action potential mode 


C 


1 sEP, 2 uEP, sLC 


Bistable mode with a stable firing and a 






quiescent state 


D 


1 uEP, sLC 


Stable firing mode 


E 


1 sEP, sLC, uLC 


Bistable mode with a stable firing and a 






quiescent state in depolarized potential 


F 


1 sEP 


Nonexcitable mode with damped oscillation 



sEP, stable EP; uEP, unstable EP; sLC, stable LC; uLC, unstable LC. 
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the change in the membrane excitability from Mode C 
to B at 8 mM [G] . To examine the role of the slow 
variables in terminating the burst, we constructed a 
bifurcation diagram by treating a single slow factor as a 
bifurcation parameter. The other slow variables were 
fixed at the values obtained at 14.6 s, just before the LP 
bifurcation. In Fig. 5 A, the bifurcation parameter is 
[ATP] . A stable LC (Fig. 5 A, blue lines) appeared over 
a higher range of [ATP] . During the burst period, [ATP] 
gradually decreased, as indicated by gray vertical lines 
sampled every 1 s. This monotonic decrease in [ATP] 
clearly promotes the burst termination through the mem- 
brane system approaching the LP bifurcation point. As 
demonstrated in our companion paper (Cha et al., 2011) , 
the terminating effect of [ATP] was mediated through 
gradual activation of outward Iratp- In Fig. 5 B, the bifur- 
cation diagram was calculated with respect to the ultra- 
slow inactivation gate of Icav (Lis) • The value of f us also 
decreased toward the LP bifurcation point, indicating 
that f us was also a key factor in the burst termination. 
Namely, the decrease of f us reduces the inward Icav and 
leads to gradual hyperpolarization. The concurrent re- 
duction of Ca 2+ influx through I Ca v lowers the plateau 
level of [Ca 2+ ];, overcoming the opposite effect of the 
accumulation of Ca 2+ in the ER. After [Ca 2+ ] , reached a 
peak at 12 s, it facilitates the burst termination (Fig. 5 C) by 
decreasing inward LjaCa and Itrpm- The accumulation of 
intracellular Na + also had a significant effect on the ter- 
mination of the burst through I NaK at 8 mM [G] (Fig. 5 D) . 



In contrast, [K + ] ; and L (inactivation state of iNaCa) 
have minor and opposite effects on the mode change 
(Fig. 5, E and F) . It should be noted that the relative im- 
portance of the slow factors on membrane excitability 
may be altered at a different [G] . 

DISCUSSION 

In the new model developed in our companion paper 
(Cha et al., 2011), the dynamics of a pancreatic (3 cell 
were described with 18 differential equations containing 
individual functional components. The present study 
focused on solving the differential equations directly, 
using bifurcation analysis. Examination was performed 
at two different levels: the entire system (Fig. 1 ) to inves- 
tigate the whole cell response to varying [G] , and at the 
level of the fast membrane subsystem (Figs. 2, 4, and 5) 
to explore the time-dependent changes in membrane 
excitability at a given [G] . The results of bifurcation 
analysis were supplemented by time-based simulations 
(Fig. 3 and Fig. 2 in Cha et al., 2011) to deduce physio- 
logically relevant conclusions. Based on the number and 
stability of the steady-state solutions (EPs and LCs) , we 
discriminated three different regions depending on [G]: 
quiescent, bursting, and continuous firing activity (Fig. 1 ) . 
We also discriminated six kinds of modes of membrane ex- 
citability depending on slow cytosolic factors (Fig. 2) . The 
exact timing of transition between modes of membrane 
excitability was also indicated on time-based simulation. 
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Figure 4. Time-dependent mode changes in 
membrane excitability during one burst-inter- 
burst cycle. (A) Time-based simulation of V m at 
8 mM [G] (gray continuous line) . Colored dots 
are EP^ EP 2 , EP 3 , and LC (min and max) in the 
membrane system, which were measured from 
the bifurcation diagrams calculated with fixed 
values of [S];S, [Ca 2+ ] ER , f us , and Ii at correspond- 
ing time points. Stable EPs, unstable EPs, stable 
LCs, and unstable LCs are indicated by black, red, 
blue, and yellow dots, respectively. The unstable 
LC (yellow) was only observed at the moment of 
mode change from Mode B to C (19.5 s) during 
<0.1 s, whereas it was not observed at the switch 
from Mode C to B during the burst (15.5 s). 
During the burst period, two sets of EPs and LCs 
were demonstrated at the sequential maximum 
or minimum of [Ca 2+ ]i. These two EP 3 s are al- 
most superimposed in the figure. The mode of 
membrane excitability at each moment is indi- 
cated at the top. (Inset) Trace of [Ca 2+ ]j during 
the same burst-interburst cycle. (B) Bifurcation 
diagrams at six representative time points in A, as 
indicated in the top left part of the figure. Dur- 
ing the burst, three time points were selected 
from those of sequential minimums of [Ca 2+ ] ; 
(11, 14.6, and 15.5 s). The same color codes were 
used for the dots as those in A. Black vertical lines 
were drawn at P Ca v = 48.9 pA mM -1 , the standard 
value in the p-cell model. 
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We investigated the ionic mechanisms for the initiation 
of the burst with lead potential analysis in our compan- 
ion paper (Cha et al., 2011), and those for the termina- 
tion of the burst with bifurcation analysis in this paper. 
Thus, the two papers together provide a novel mathemati- 
cal description about the (3-cell bursting activity and the 
role of individual functional units or molecules in the 
response to glucose. 

Comparison with previous studies in respect 
to bifurcation analysis 

In the past few decades, bifurcation analysis has been 
used to clarify the mathematical principles underlying 
bursting activity in pancreatic (3 cells. In most studies, 
extremely simplified models have been used because 
they were amenable to mathematical analysis. In this 
study, we applied bifurcation analysis to a complex 
(3-cell model based on extensive experimental data, 
with the aim of clarifying important physiological 
mechanisms in reference to individual functional com- 
ponents. We demonstrated that the bursting activity in 
this complex system is generated with the same basic 
principles as before. 

Interactions between fast and slow systems. Previous mod- 
eling studies consistently indicated that the bursting 
rhythm is generated by the reciprocal interactions be- 
tween fast and slow subsystems. These studies presented 



useful hypotheses for a negative feedback mechanism 
between an ionic current with one single slow variable. 
For example, early (3-cell models hypothesized an inter- 
action between [Ca 2+ ]i and Ca 2+ -dependent activation of 
Ikce (Chay and Keizer, 1983; Sherman et al., 1988) or 
Ca 2+ -dependent inactivation of Icav (ChayandKang, 1988). 
Subsequendy, other hypothetical mechanisms include 
voltage-dependent slow inactivation of Ic a (Chay, 1990; 
Keizer and Smolen, 1991), [ATP], or [ADP] to modu- 
late the conductance of Ikatp (Smolen and Keizer, 1992) , 
or [Ca 2+ ] ER to activate I SO c (Chay, 1996, 1997). With our 
detailed cell model, we demonstrated that the following 
multiple slow variables work in concert to generate the 
burst-interburst rhythm: slow inactivation (f us ) via Icav* 
intracellular ATP metabolism via Ikatp, [Ca 2+ ]j via I^ca 
and Itrpm! and [Na + ]; via I^aK- These mechanisms make 
a positive contribution to termination of the burst against 
the minor and opposite effects of [K + ]j or Ii (inactivation 
state of INaCa) • Moreover, our model showed that even a 
single slow variable has complex influences on the mem- 
brane excitability. For example, [ATP] promotes the 
burst termination by activating Ikatp but, at the same 
time, retards it by depressing I Na K; [Ca 2+ ]i stabilizes the 
oscillation during the initial phase of the burst but ter- 
minates the burst in the late phase. 

Transitions between an EP and an LC. Results in this study 
are consistent with the conclusions of previous studies, 
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Figure 5. Effects of individual slow variables on the mode change toward termination of the burst. Each bifurcation diagram was ob- 
tained by varying single slow variable, [ATP] (A),f us (B), [Ca 2+ ]i (C), [Na + ]; (D), [K + ] ; (E),andli (F) as a bifurcation parameter, with the 
remaining slow variables fixed at their values at 0.9 s before the LP bifurcation (open circles) . Stable EPs, unstable EPs, stable LCs, and 
unstable LCs are indicated by black, red, blue, and yellow lines, respectively. Six gray vertical lines in A-D indicate the values of the cor- 
responding slow variables sampled at an interval of 1 s (from 11 to 16 s in Fig. 4 A) . In C, the sequential values of the minimum [Ca 2+ ] j 
were sampled. In E and F, two gray lines indicate the corresponding values at 1 1 and 16 s. Arrows represent the sampling sequence. 
We confirmed that the bifurcation diagrams were qualitatively the same when slow variables were fixed at different time points during 
the burst. 
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namely that burst-interburst rhythm at a given [G] can 
be explained by repetitive transitions between a stable 
EP and a stable LC, although our graphical analysis is 
largely different. In the most successful and straightfor- 
ward presentation using a simple model (Sherman et al., 
1988) , a bifurcation diagram was obtained using [Ca 2+ ]; 
as the bifurcation parameter, which was the sole slow 
variable in their model. The bifurcation diagram was 
superimposed with a Ca 2+ nullcline in V m -[Ca 2+ ]; space, 
and thereby, the overall behavior of the system was eas- 
ily tracked along the V m -[Ca 2+ ] ; diagram guided by both 
the Ca 2+ nullcline and bifurcation points. In the study of 
Bertram and Sherman (2004), their model was com- 
posed of three slow variables, [Ca 2+ ] ; , [Ca 2+ ] E R, and [ADP], 
and bifurcation diagrams were presented with respect 
to [Ca 2+ ];, by fixing the other two variables. Therefore, 
the combined effects of multiple slow variables were 
only indirectly inferred by comparing bifurcation dia- 
grams calculated with two different values of [ADP] or 
[Ca 2+ ]£R- Our model is an even more complex system, 
with eight slow variables, so it was extremely difficult to 
prepare such a straightforward bifurcation diagram. To 
overcome this problem, we developed an alternative way 
to show the transition of V m between an EP and an LC 
along the time axis. To show the net effects of the slow 
variables during the bursting rhythm, we used the val- 
ues of all the slow variables at each time point to calcu- 
late EPs and LCs. In our presentation, time-dependent 
changes in the mode of membrane excitability were ex- 
plicitly identified on the record of the time-based simu- 
lation (Fig. 4 A). 

Modal changes in behavior of the whole system. We found 
that the bursting activity ceased by decreasing [G], and 
the firing became uninterrupted by increasing [G] in 
our (3-cell model. These transitions from quiescence to 
bursting or to continuous firing in the whole system 
were consistently observed in several previous studies, 
but only indirectiy. Namely, changes in glucose were 
mimicked by increasing either the rate of cytosolic Ca 2+ 
sequestration (Chay and Keizer, 1983) or the Ca 2+ - 
binding affinity to Ca 2+ channels (Chay, 1993). In more 
complex models, a hypothetical parameter proportional 
to the proton motive force (Keizer and Magnus, 1989) 
or the conductance of Iratp (Smolen and Keizer, 1992; 
Bertram et al., 1995) was changed, instead of calcula- 
ting the reaction pathways for the glucose signal. There- 
fore, these simple models did not directly address the 
key question of how changes in [G] modulate cellular 
activity. In this study, we simulated the whole spectrum 
of [G] dependency by implementing the changes in 
metabolic status after changes in [G] (Fig. 1). This en- 
abled us to estimate the values of [G] at which the 
cell changes its electrophysiological characteristics. 
Although we successfully reproduced cellular response 
to a range of [G] relevant to experimental measurements 



in the mouse, improvements in the metabolic compo- 
nents are still awaited to get a deeper insight into the ef- 
fects of glucose. In the future, we would especially like 
to examine the effects of [Ca 2+ ]i on the production of 
ATP through the tricarboxylic acid cycle and oxidative 
phosphorylation. 

A wide range of cycle length in bursting activity. In a simple 
model possessing a single slow variable, the range of 
burst cycle length was rather limited if compared with 
the experimental observations. Bertram et al. (2000) 
demonstrated that burst cycle length can be varied over 
a wider range by including two different slow processes, 
one with a relatively small time constant (si) and the 
other with a much larger time constant (S2) . Three slow 
variables were assigned to [Ca 2+ ];and [ATP] or [Ca 2+ ]ER 
in a subsequent study (Bertram and Sherman, 2004) . In 
the present study, we used nine slow variables, includ- 
ing seven substrate concentrations as well as the slow 
gating of Icav and iNaCa- This resulted in a wide range of 
burst cycle length when external [G] was changed (Fig. 2 
in Cha et al., 2011). Based on the theory of Bertram etal. 
(2000), comparing the time courses of the slow vari- 
ables implies that [ATP] ([MgADP]) or f us might corre- 
spond to Si, and [Na + ]; or [Ca 2+ ]£R to S2, in our model. 
That is, the duration of one cycle of the bursting 
rhythm is short at a low [G] (8 mM) where the varia- 
tions in [ATP] or f us were predominant, whereas rela- 
tively large variations in [Na + ]; or [Ca 2+ ]ER govern a 
much slower rhythm at a high [G] (16 mM). It should 
be noted that the rate of change in [ATP] was largely 
affected by the Ca 2+ -dependent consumption of ATP 
(k AT p,ca; Eq. S103 in Cha et al., 201 1) , and that of [Na + ] ; 
was determined by NaK or NaCa activity in our model. 
Thus, for more reliable reproduction of the experi- 
mental burst duration, the above parameters should be 
improved in the future, when more extensive experi- 
mental data are available. 

Conclusion 

In conclusion, collectively with the time-based integra- 
tions, steady-state solutions of our differential equations 
explicitly proved the hypothesis about the roles of indi- 
vidual ion channels and transporters in generating 
(3-cell electrical activity in relation to their complex inter- 
actions with slow variables. Furthermore, working hy- 
potheses for new experiments can be obtained from a 
mathematical model with detailed membrane compo- 
nents and cytosolic mechanisms. Although quantitative 
predictions from any mathematical model are depen- 
dent on how correctly the individual model components 
are described, this study and our companion paper (Cha 
et al., 2011), using bifurcation analysis, lead potential 
analysis, and time-based simulations, provide a frame- 
work to improve an objective understanding of this 
complex system. 
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APPENDIX 
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By combining Eqs. A4 and A5, 
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By combining Eqs. Al, A2, A3, and A6, 



d [Ca 2+ ] ER vol f ER ( dV„ _ dfNa^ _ d[K+f _ 2 d[Ca 2 *] t ) 
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By integrating both terms with t from t = 0, 



9 . vol. f CD C 
[Ca 2+ ] ER = —^5. { _^L_ (Vm - V n (0) ) - ( [Na + ], - [Na + ], (0) ) ( A8) 

2 vol EE volj F 

-( [K + J, - [K + ], (0) ) - - ( [Ca 2+ ], - [Ca 2+ ],(0) )} + [Ca 2+ ] ER (0) . 
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